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ABSTRACT 

Two dimensional particle disks in proximity to a planet are numerically integrated to 
determine when a planet in a circular orbit can truncate a particle disk. Collisions are 
treated by giving each particle a series of velocity perturbations during the integration. 
We estimate the mass of a planet required to truncate a particle disk as a function of 
collision rate, related to the disk optical depth, and velocity perturbation size, related 
to the disk velocity dispersion. We find that for particle disks in the regime estimated 
for debris disks, a Neptune mass planet is sufficiently massive to truncate the disk. If 
both the velocity dispersion and the disk optical depth are low (dispersion less than 
approximately 0.02 in units of circular motion, and optical depth less than 10 -4 ) then 
an Earth mass planet suffices. We find that the disk is smooth and axisymmetric unless 
the velocity perturbation is small and the planet mass is of order or greater than a 
Neptune mass in which case azimuthal structure is seen near prominent mean motion 
resonances. 



1 INTRODUCTION 



Recent observations made with the Spitzer Space Tele- 
scope have added to the total number of known debris 



and circumstellar disks (lRiekell2005l ; iBeichman et al.l 



Gorlova et all 1 20061; IChen et alj 120061 ; iBrvden et ail 



2005; 



2006; 



Beichman et al. 20061) . Theses disks have total infrared 



the disk, though when possible, direct imagi ng in the in- 
frare d can provide better measurements (e.g.. iMarsh et al.l 
120051 ). 

Planets have been proposed to explain t he morphology 
of a number of systems with dusty disks (e.g..lOzernoy et al. 
20001; IWilner et al.112002; Qui llen fc Thorndikell2002l ; Iwvatt 



opacity, as estimated from the fraction of stellar light re- 
emitted in the infrared, in the range tir ~ 10~ 3 — 10 -5 . 
Infrared spectra of many debris are well fit with a single 
black body temperatur e suggesting that they possess inner 
holes (|Chen et al.ll2006T ). Planets are suspected to be respon- 
sible for the inner clea rings (e.g. jKalas et al.ll2005l ; iQuillenl 
l2006l ; IChen et al.ll2006l ). 

The fraction of starlight absorbed by the dust is re- 
lated to the timescale between dust particle collisions. Here 
we focus on collisions as they are expected to make signif- 
icant modifications to the dynamics. The relation between 
the fraction of starlight absorbed by the dust and the colli- 
sion timescale depends on the width and radius of the disk. 
A ring with radius r, width dr, and optical depth normal 
to the disk plane, r„, has total surface area A — 2nrdrr n 
however at radius r the starlight passes through a sphere of 
area A-kt 2 . For low opacity systems the infrared opacity is 
related to the disk optical depth by 



20031 ; iDeller fc M addison 200jJ), however most models do 



not take into account collisions between particles in the disk. 
The collision timescale, t co i, for the smaller particles can be 
estimated from the optical depth, r n , normal to the disk 
plane, 



t c 



(3r„n) 



2r 
dr 



tir. 



(1) 



We expect that the disk optical depth, r n , exceeds tir by 
a modest factor (dr < r) and that r n is the same order of 
magnitude as tir unless the ring is very narrow. The in- 
frared opacity measurements allow us to crudely estimate 
the normal disk optical depth and so collision timescale in 



(2) 

IIHanninen fc Sa iol ll992D where n is the mean motion at ra- 
dius r . As emphasized bv lDominik. fc Decinl l|2003l ); IWvattl 

(2005) , for most debris disks the collision timescale be- 
tween particles is shorter than the Poynting-Robertson drag 
timescale and so collisions should be considered wh en inter- 
pretin g the evolution and morphology of these disks. iQuillenl 

(2006) proposed that the steep inner edge of Fomalhaut's 
disk might be due to the removal (by the planet) of parti- 
cles sca t tered interior to the disk edge by collisions. While 
Quillen (2006) suggested that the chaotic zone near the coro- 
tation region is a likely location for disk truncation, this has 
yet to be tested by numerical simulations that incorporate 
both perturbations by a planet and collisions on a timescale 
consistent with that expected for debris disks. That is what 
we aim to do here. 

For disk optical depths in the regime of debris disks, 
T n ^ 10 -2 , spiral density waves cannot be driven at Lind- 
blad resonances (|Quillenl [2006; this regi me is also dis- 
cussed in the context o f plan e tary rings by iFranklin et al.l 
Il98d ; lHanninen fc Said Il992l ; iLissauer fc Espresatd Il998; 
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Figure la. Azimuthally averaged density as a function of ra- 
dius for simulations C1,D1,E1,F1,M1,N1 and LI listed in table 
[T] These simulations have different numbers of collisions per or- 
bit and velocity perturbations, but the same planet mass ratio, 
fj, = 10 -4 . Radii are given in units of the planet's semi-major 
axis. Arrows are drawn at the radii corresponding to the location 
of the chaos zone boundary or at r = 1 + 1.5/x 2 / 7 for each planet 
mass ratio that has a displayed simulation. 



lEspresate fc Lissauerll200ll ). However collision induced an- 
gular momentum transfer causes a thin disk to diffuse radi- 
ally and be come wider. A particle ring has two characteristic 
timescales (|Brahidll977h . the collision timescale, t co i, that is 
related to the disk optical depth, and the diffusion timescale, 
tdiff, that depends on both the disk velocity dispersion and 
collision timescale. 

The radial evolution of an isolated particle ring can be 
described with a diffusion equation 



dN__ d_ ( D dN_\ 



(3) 



where N(r) is the number d e nsity as a function of 
radius r, (e.g.. | Petit fc Henonl Il987l equation 48 and 



iLithwick fc Chiand l200rj equation 11). The diffusion coef- 
ficient, D, depends on the collision timescale and velocity 
dispersion, u 2 , in the disk, 



(4) 

This diffusion coefficient is similar to a viscosity and can 
be estimated by considering the mean free path and parti- 
cle velocity differences se t by the epicyclic amplitude (e.g., 
iMurrav fc Dermottlll999T ). As the collision time depends on 
the number density, the diffusion coefficient can d epend on 
N (|Petit fc Henonll 19871 ; ITithwick fc Chiand |2006h . 

In this paper we adopt a diffusive approach toward sim- 
ulating the affect of collisions on a particle disk that re- 
sides near a planet. When the disk optical depth is low, 
many orbits pass between collisions. During this time par- 
ticles can be integrated directly under the influence of 
only gravity. Previous works simulating collisions adopt a 
statistical method for predicting wh en each particle suf- 
fers a collision and its affect (e.g., lEspresate fc Lissauerl 
l200ll ; iMelita fc Woolfsonl Il998l ). or search for close ap- 
proaches between particles and t hen compute momen- 
tum changes for both particles (e.g. iLewis fc Stewardl200fi 



Figure lb. Similar to Hal except simulations C2,D2,E2,F2,M2, 
and N2 are shown with a planet mass ratio of fi = 10 -5 . 
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Figure lc. Similar to Hal except simulations C3,D3,E3,F3,M3, 
and N3 are shown with a planet mass ratio of fj, = 10 -6 . 



ICharnoz et al]|200ll ; ILithwick fc Chiandl2006h . In this work 
we adopt the first approach. We only perturb particles 
on a mean collision timescale and compute a pertur- 
bation to the particle's momentum. Both the collision 
timescale and momentum perturbation are chosen indepen- 
de nt of the particle distributi on (similar to the approach 
by lEspresate fc Liss aucr 2001). A more sophisticated sim- 
ulation would allow the particle distrib ution to set both 
quantities (e.g.. IMelita fc Woo lfson 1998]). However this re- 
quires additional computational effort, either by com put- 
ing close approaches carefully l|Charnoz et al.l [200lj) , let- 
ting the computed collisi ons depend on a numerical grid 
ILithwick fc Ch iang 2006) or integrating sufficient numbers 
of particles that the velo city distribution is well sampled 
jMehta fc Woolfsonl Il998r i. 



2 NUMERICAL INTEGRATIONS 

Numerical integrations were carried out in the plane, using 
massless particles under the gravitational influence of only 



the star and a massive planet with zero eccentricity, e p 



0. 
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and semi- major axis, a p = 1, using a conventional Burlisch- 
Stoer numerical scheme. Initial particle mean anomalies and 
longitudes of perihelia were randomly chosen. In each sim- 
ulation 1000 particles were integrated for ~ 10 5 planetary 
orbital periods. 

Each particle suffers collisions with a mean collision rate 
set by the collision timescale. The parameter we use to de- 
scribe the collision rate is the number of collisions per planet 
orbit, N c , for an individual particle; N c w P/t co i, where P is 
an planet orbital period. The parameter N c is related to the 
norm al disk optical depth by N c ~ 18r n (|Hanninen fe Salol 
1992). Each planetary orbit the number of collisions is com- 
puted from the total number of particles in the simulation 
and the collision rate. The particles that receive collisions 
are chosen randomly from those integrated. These particles 
receive perturbations to their momentum at a randomly cho- 
sen time during the orbit. 

If a particle suffers a collision its radial velocity is ad- 
justed to damp the energy, 

v r — » ev r (5) 

where —1 < e ^ 1. Here we have set e = 0. Had we not re- 
stricted our simulations to zero eccentricity planets we would 
have required the radial velocity perturbation to depend on 
the mean set by the distribution of particles. The tangential 
velocity is given a random kick 



V8, 



VK 



(7) 



ity — > ity + Av 



(6) 



where Av is chosen from a normal distribution with disper- 
sion Odv 

Because the velocity perturbation at each collision, Av, 
is chosen from a normal distribution its expectation is zero. 
While total angular momentum should not drift rapidly, an- 
gular momentum is not conserved during the simulation. 
Adjustment or redistribution of angular momentum follow- 
ing collisions would be required to co nserve angular mo- 
ment um (e.g., as explicitly described bv lMelita fc Woolfsonl 
1998). We do not chose to do this as we continually regen- 
erate particles that have scattered outside the region of in- 
terest. To conserve angular momentum computational time 
would have been wasted integrating particles that diffuse 
away from the planet as well as toward it. Our adopted pro- 
cedure allows us to simulate a diffusive disk with minimal 
time spent integrating particles outside the region near the 
planet. 

Particle initial semi-major axes were chosen to lie be- 
tween 1.4 and 1.5 times the semi-major axis of the planet. 
This way the particles were allowed to diffuse from larger 
radius, approaching the planet. Particles are removed from 
the simulation if they have a semi-major axis greater than 
1.6 or less than 0.7, have an eccentricity greater than 0.5 or 
are unbound and if they move closer to the planet than 1/10 
of its Hill radius. Particles that are removed from the sim- 
ulation are regenerated choosing their orbital parameters in 
the same way as the initial particles were chosen. This en- 
sures that a constant number of particles is integrated at all 
times. 

The initial particle eccentricities were drawn from a 
Rayleigh distribution with a mean eccentricity of emit that 
is chosen to match the velocity perturbation parameter, adv 
or 



where vk is the velocity of a particle undergoing circular 
motion at the radius of the planet. A particle in a nearly cir- 
cular orbit with eccentricity e has radial velocity component 
= e sin nt and tangential component = ( 1 + 1 cos ni) . 
A distribution of particles with the same eccentricity and 
random mean longitudes and longitudes of periastron has 

<vl> _ £_ 

8 ■ 



radial dispersion < "j ,> 



and <W-*> 2 > 



2 „2 iv z 

It. K K 

This last expression accounts for the factor of Vo in equa- 
tion[7] We describe and refer to simulations in terms of either 
the velocity perturbation, adv, or the mean initial particle 
eccentricity, emit, as they are directly related to one another. 

A steady state is reached after the disk has had time 
diffuse from the generation region (at a radius of 1.4-1.5) 
to the particle removal region, near the planet's orbital ra- 
dius. We checked that the radial velocity dispersion distant 
from the planet (greater than r ~ 1.2) is that expected 
from perturbations given to the tangential velocity pertur- 
bation at each collision or < v 2 >f» 4<r^. This is consistent 
with a distribution of particles with a mean eccentricity of 
< e >= emit = \'' iv ■ Because the radial velocity compo- 
nent is damped at each collision the mean particle eccen- 
tricity does not significantly increase past this value except 
near the planet. 

A density distribution is created from particle positions 
at the same time (planet mean anomaly) in each planetary 
orbit. We have added 1000 of these from different timesteps 
(each at the same time in the planet's orbit) to ensure that 
the density distributions are smooth and have a large dy- 
namic range. We work in units of the planet's semi-major 
axis, a p , and orbital period. The mass of the planet is de- 
scribed in terms of its mass ratio, \x, the ratio of the planet 
mass to that of the central star. We describe the velocity 
perturbation aj v m units of the velocity of a particle in a 
circular orbit, vk, or in terms of the associated initial par- 
ticle eccentricity, emit- 



3 SIMULATIONS 

The simulations we discuss have parameters listed in Table 
[1] We relate the simulation parameters to quantities that 
can be constrained by observations, the disk normal optical 
depth, r n , and the velocity dispersion that is related to the 
disk thickness. We have run simulations for three different 
planet mass ratios, \i — 10 -4 , 10 -5 and 10 -6 , two different 
collision rates, N c — 10 -2 and 10 -3 collisions per parti- 
cle per orbit corresponding to disk normal optical depths 
t„ = 0.5 x 10~ 3 and 0.5 x 10 -4 and four different veloc- 
ity perturbation sizes. In Table [1] the velocity perturbation 
sizes are listed in terms of the initial mean eccentricity ei n u 
where (Jdv/vK = einit/V8 and are related to the disk veloc- 
ity dispersion by {u/v K ) 2 ~ e 2 lnit /2. 

Azimuthally averaged density profiles (number density 
as a function of radius) after equilibrium is reached are 
shown in Figures [la] - [4a] where each figure shows simula- 
tions with only one planet mass. A comparison of Figure [Tal 
- [Tc] suggests that the radius at which the particle number 
density drops is related to the planet mass. 

In the collisionless restricted 3 body system there is 
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Figure 2a. Number density as a function of r and <f> for simula- 
tion C2 with planet mass ratio fi = 10 -5 . Little structure is seen 
in the disk. 

an abrupt change in dynamics in the corotation region at 
a semi-major axis that we refer to as the boundary of 
the chaotic zone. The width of this zone has been mea- 
sured numerically and predicted theoretically by predict- 
ing the semi-major axis at which the first order mean mo- 
tion resonances overlap (IWisdomlll980l ; iDuncan et al. I ll989l ; 
iMurrav fc Holmanl [l997h . In our figures we have with an 
arrow marked the chaotic zone boundary at a radius r z = 
a p + da z with the dista nce between the planet and chaotic 
zone da z = 1.5a p /i 2/7 \ Wisdom! liliol ; IDuncan et al.lll989T l 
for the different simulations. We see from the figures that 
this boundary is approximately the location at which the 
radial profiles suffer a drop near the planet, as would be ex- 
pected from the change i n the dynamics due to overlap of 
mean motion resonances jWisdomlll9"80l ; lMurrav fc Holmanl 
ll997l ; lQuillen fc Faberll2006l ). 

We note that the density drops before the chaotic zone 
boundary for the simulations with the most massive planets. 
The simulations with the most massive planet and the steep- 
est profiles are also the ones that show the most structure. 
Figures l2"al and [2bl display density distributions as a function 
of radius and angle for simulations C2 and Fl. We find for 
planet mass ratios equal to or below 10 -5 little structure 
is evident in the density distribution. However for a planet 
mass ratio of 10 -4 and low values of the velocity perturba- 
tion, structure at strong mean motion resonances such as the 
3:2 and 3:4 resonance is present. When the resonances are 
strong trajectories can become planet orbit crossing. Par- 
ticles that diffuse into these regions are removed from the 
simulation causing structure in the surface density (as seen 
in figure [2b)l and a reduction in the density profile at a ra- 
dius larger than the chaotic zone boundary. This is likely to 
account for the density drop past the chaotic zone boundary 
evident in the profiles shown in figure [Tal 

In Figure [3] we show profiles for the 3 simulations with 
the same collision rate and perturbation but different planet 
masses and with radius rescaled by ii 2 ^ 7 so that the [i = 10~ 4 
and fj, = 10 -6 simulation chaotic zone boundaries lie on 
that for the fj, — 10 -5 simulation. In other words r has been 
rescaled in this figure by a factor of (10 _5 //i) 2//7 . Here we see 
that this radius does give an estimate for the radius at which 
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Figure 2b. Number density as a function of r and <f> for simu- 
lation Fl with planet mass ratio [i = 10 -3 . Structure is seen in 
this disk associated with mean motion resonances. The 2:3 mean 
motion resonance is located at r = 1.31 and the 3:4 mean motion 
resonances at 1.21. Both resonances are associated with structure 
in the disk. 

a disk is likely to suffer a drop in density due to removal of 
particles by the planet. The rescaled profiles do not lie on 
top of each other; the lower mass planet simulations have 
shallower profiles after rescaling. The disk slopes are not self 
similar after rescaling implying that the disk profile shape 
is not solely a function of the distance to the chaotic zone. 

The slope of the disk edge is dependent on the collision 
rate and perturbation size as well as the planet mass. In 
figures I4al I4bl and |4cJ we show profiles for the same num- 
ber of collisions per orbit and collision velocity perturba- 
tion sizes but for different planet masses. The simulations 
shown in I4bl have twice the velocity perturbation size and 
the same collision rate as those shown in 2a] whereas those 
shown in [4c] have the same velocity perturbation size but 
a tenth the number of collisions per orbit as those shown 
in I4al We find that a larger change in the profile shape is 
caused by variation in the velocity perturbation parameter 
<Jdv than variation by the same factor in the collision rate 
No. 

3.1 Diffusive descriptions for the density drop 

We now discuss diffusive descriptions for the disk edge to 
better understand the parameters affecting the decrease in 
particle density near the planet. Since the velocity dispersion 
is u 2 ~< v 2 > and < v 2 >~ 4<r 2 .„ we estimate that the 
diffusion coefficient (equation [3} is 



In units of the square of the planet semi-major axis divided 
by the rotation period, 

D ~mA—\ 2 = ^i^. (9) 

V VK J 2 

The diffusion equation for an isolated ring, Equation [31 
must be modified to include the effect of the planet. When 
a steady state is reached (^ « 0), a simple model is 
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Figure 3. Azimuthally averaged density as a function of radius 
for simulations with different planet mass ratios C1,C2, and C3 
but the same collision rate and velocity perturbation size. Here 
the radius has been rescaled by a factor of (10 — 5 //x) 2/7 so the 
chaotic zone boundaries for the three simulations line up. 
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Figure 4a. Azimuthally averaged density as a function of radius 
for simulations with different planet mass ratios but the same 
collision rate and velocity perturbation sizes. Here simulations 
C1,C2, and C3 are shown. 



(10) 



d_ ( D 9N_\ _ dN_ 

dr \ dr J dt 

where the right hand side represents the 
rate that particles are removed from radius r because of 
perturbations from the planet. This term is expected to be 
a function of planet mass and distance to the planet but 
could also be a function of the particle velocity distribution 
or eccentricity distribution. 

We assume that there is only one important scale in the 
problem, the distance between the planet and the chaotic 
zone boundary, da z = 1.5/i 2//7 . We rewrite equation 1 101 with 
in terms of variable y — - " d ^ v and assuming that the diffu- 
sion coefficient is independent of r and N, 



2 N 



dy 2 



D 



dN 
~dt 



(11) 



A simple form for the right hand side of equation 1101 de- 
scribing removal of particles by the planet, would be 



dN 
~dt 



planet 



Nf(y) 







for < y < 1 
for y ^ 1 



When f(jj) = 1 the solution 
N{y) oc e yl 

with inverse scale length 



/ ■ 



(12) 



(13) 



(14) 



11 f(y) = 1 — 3/, then N(y) is an Airy function and if 
f(y) = exp(— y) (as might be expected from lifetime mea - 
surements in collisionless integrations; iDavid et al 1 12003T ) 
the solution is a modified Bessel function. These solutions 
also behave in an exponential manner with an inverse scale 
length that is of the same order of magnitude as /. Figures 
Hal - [Tc] [4a] - [4c] show that the density profiles do approach 
an exponential form near the planet where the density gra- 
dient is steepest. In parameters related to our simulations 
we expect the inverse scale length 
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Figure 4b. Same as figure l4al but for simulations N1,N2,N3. 
These simulations have twice the velocity perturbation size and 
the same collision rate as those shown in !4al 



(15) 



From the simulations we measure the density near the 
planet compared to that away from the disk edge or specif- 
ically the log of this ratio 



log 



( n(v = Q) 

\N(y>l) 



(16) 



Equations 1131 and 1151 suggest that A depends on parameters 
describing our simulations, such as the planet mass ratio, 
collision rate and velocity perturbation size, however the 
above approximations and assumptions may be too simplis- 
tic to give an accurate prediction for the inverse scale length 
I. Furthermore I contains the unknown timescale, t rernove . 

The density decrement, A, measured from the radial 
profiles, is plotted in Figures [5a] and I5bl As predicted by 
equation 1151 our measured density decrement is consistent 
with a power law function of the simulation parameters. Also 
as we expected A increases with increasing planet mass ra- 
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Figure 4c. Same as figure l4al but for simulations D1,D2,D3. 
These simulations have the same velocity perturbation size but a 
tenth the number of collisions per orbit as those shown in I4al 



tio, and decreases with increasing collisions rate and veloc- 
ity perturbation size. To measure the exponents, we have fit 
lines to the density decrements finding pretty good corre- 
spondence with the following function that is also plotted in 
Figures [5a] and I5bl 



logio A 



0.12 + 0. 23 log 10 
-0.451og 10 



(^-""-(sr) 

VaoT/ 



io- 2 

(17, 



This function implies that the inverse scale length 
/ u \ ' 23 / N \"°- 1 /p- t \-°' 45 

'~ 3 -°(ifO GfO im) 



(18) 



A comparison between the measured scale length above 
(equation [T8]) and the predicted by equation [15] reveals some 
inconsistencies. Equation [15] implies that the length scale 
should be proportional to Na X e^a rather than N~°' 1 e'init 
as we find above. The above function measured I combined 
with equation [15] would imply that 



remove 



145 ! 



(&r (M mr <»> 



in orbital periods. It is not clear why the removal timescale 
would increase with planet mass. Had we used a larger scale 
for the distance than the distance to the chaotic zone bound- 
ary in equations [T2] and [T3] then the form of t remove would 
have changed. It is expected that tremove depends on the 
velocity perturbation size as we do expect the particle life- 
time near the planet to depend on particle eccentricity, with 
shorter lifetimes or removal timescales at higher particles 
eccentricities. Perhaps the particle lifetime depends on the 
collision rate because particles tend to diffuse to higher and 
higher eccentricity between collisions. Our measured I is not 
as strong a function of N c as predicted by equation 1151 A 
better theoretical model for both particle lifetime and diffu- 
sion is required to better match the measured exponents of 
equation [18] 



3.2 Disk truncation 

The shallowest profiles displayed in Figures [la] - llbl repre- 
sent particle disks that have sufficiently high collision rates 
and velocity dispersions that a planet cannot efficiently trun- 
cate the disk. It is interesting to apply our measured density 
decrement to estimate the minimum planet mass that can 
truncate a disk. We set a limit of log 10 A = 0.12 (correspond- 
ing to a density drop of 20) from the shallowest simulations 
and solve equation [T7] to find the mass ratio of a planet that 
cannot effectively truncate at disk of a particular collision 
rate and velocity dispersion. To truncate a low optical depth 
particle disk a planet must have mass ratio greater than 

log 10 I* > -6 + 0.43 log 10 (i) + 1.95 log 10 (g^) (20) 
or in quantities more directly related to observations 



log 10 M> -6+0.43 log 10 ( 



5 x 10" 3 



The function (equation I20[) is plotted in figure [6] for col- 
lision rates N c — 10~ 2 and 10 -3 collisions per orbit cor- 
responding to disk optical depths t„ — 0.5 x 10~ 3 and 
0.5 x 10~ 4 . Instead of plotting this minimum mass ratio as a 
function of dnit we plot as a function of velocity dispersion 
using ef nit ~ 2(u/vk) 2 - On this plot we also show a limit 
corresponding to the eccentricity of a particle that would 
allow a particle to cross from the chaos zone boundary to 
the planet's semi-major axis, or the limit (u/vk) ~ 1.5/j, 2 ^ 7 . 
We refer to this limit on the plot as /x e - 

It is interesting to compare our mass limit for disk trun- 
cation to that expected for an a accretion disk. A gaps opens 
in a viscous accretion disk when the torque driven by spiral 
density waves excited by a planet exceeds viscous inflow or 

H > 407?e _1 (22) 

(|Lin fc Papaloizoul Il979l ; iBrvden et afl Il999l ; IWardl Il997l ). 
Here the inverse of the Reynold's number, Re^ 1 = 
a(u/vi() 2 - We have plotted this function on Figure [6] for 
a — 0.001. A particle disk has effective viscosity set by the 
diffus ion coefficient of equation |41 (e.g.. iMurrav fc Dermottl 
1999), hence if our particle disks behaved similar to an ac- 
cretion disk then this line (shown as a thin dotted line in 
figure [SJ would like approximately on top of the dotted line 
for 7V C = 10 -2 with a ~ N c /6. It is interesting that the 
disk truncation criterion scales with velocity dispersion in 
the same way as the accretion disk; both require that the 
planet mass ratio exceed a constant times the square of the 
velocity dispersion. We find that a larger planet is required 
to open a gap in an accretion disk with equivalent viscosity 
to the diffusion coefficient in a low optical depth particle 
disk. This is somewhat counter-intuitive as spiral density 
waves excited by a planet carry angular momentum pushing 
material away from the planet and they cannot be driven in 
a low optical depth particle disk. However collisions damp 
orbital eccentricity and high eccentricity orbits can be more 
quickly scattered by a planet. A lower mass planet might be 
able to truncate a particle disk because the long timescale 
between collisions allows particle eccentricities to grow near 
the planet. 

Figure [5] illustrates that for particle disks in the regime 
expected for debris disks, a Neptune mass planet is suffi- 
ciently massive to truncate the disk, and if the velocity dis- 
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persion is low (u/vk < 0.05) and the disk optical depth low, 
T n < 10~ 4 , an Earth mass planet suffices. 



4 DISCUSSION 

Here we have n ot considered models with dust produce d in 
resonance fe.g jQuirlen fc Thornd"ike1l2002l ; IWvattll2003l ') but 
only allowed particles to diffuse inward toward resonances. 
In resonance, certain angles with respect to the planet have 
longer removal timescales than others. Effectively | ; t 
depends on the orbital elements of the particle distribution 
and the diffusion equation should be expanded to depend not 
on radius but all the orbital elements. Individual resonances 
may play an important role setting the velocity dispersion 
disks near larger pl anets causing increased diffusion nea r 
the disk edge (e.g., iQuillenl 12009 : IQuillen fc Fabeij I2006T ). 
We expect additional complexity and structure if particles 
are produced trapped in resonance or if the system is not in 
a steady state. The simulations were run an average of 10 5 
orbits which is long compared to the age of some large de- 
bris disks. Our integration times are long because we allowed 
particles to diffuse toward the planet from a moderately dis- 
tant location and we ran each simulation sufficiently long to 
achieve a steady state. Additional phenomena may be dis- 
covered in evolving systems. 

The simulations were carried out with a collision rate 
and associated collision induced velocity perturbation inde- 
pendent of position, however particles near the planet should 
have suffered fewer collisions because of the steep drop in 
density. Particles in the disk edge should have also suffered 
larger velocity perturbations as the velocity dispersion in- 
creases nearing the planet. The velocity dispersion increase 
(spanning a factor of a few for the largest planet) is not as 
significant as the density decrease nearing the planet (span- 
ning orders of magnitude for the largest planet simulations) . 
A more accurate simulation probably would have a steeper 
disk edge due to the decrease in collision rate near the planet 
and would have required longer to achieve equilibrium as a 
consequence. 

Here we have adopted a diffusive approach toward sim- 
ulating particle disks and neglected the role of the particle 
size distribution. If collisions are destructive then diffusion 
caused by collisions can still take place. However each diffu- 
sive step involves a destructive collision producing smaller 
bodies and these in turn have a shorter collision timescale. 
We imagine an effective diffusion coefficient that depends 
not only on local particle density but on the size distribu- 
tion with the larger particles effectively having longer colli- 
sions times and smaller diffusion coefficients. The different 
sized particles also may have different velocity dispersions. 
In this case we would expect that the larger bodies would 
have smaller diffusion coefficients as the large bodies would 
have lower velocity dispersions. It is tempting though pre- 
mature to predict that the different size particles would ex- 
hibit different disk morphologies as a consequence of their 
different collisional timescales and velocity dispersions. 

It is interesting to compare our simulations to imag- 
ing studies of debris disks. Our lower mass limit (equation 
1211 see circle in Figure [6} is n ear but a few times lower 
than the previous l imit set by IQuillenl (l2006f ) for Foma- 
lhaut. As noted by IQuillen fc Faber 1 200a ) the dynamics 



of low free eccentricity particles near an eccentric planet 
is similar to low eccentricity particles near a planet in a 
circular orbit, so we can extend equation [5U to low eccen- 
tricity planets. Observations of Fomalhaut's disk have not 
revealed any lu mps or asymmetri es in the disk, other than 
its eccentricity l|Kalas et al.l l2005l. consequently the planet 
mass must not be high (fx < 10 -4 ). Otherwise Fomalhaut's 
disk would look like the disk displayed in Figure I2bl Con- 
sequently we also confirm, in a differe nt way, the u pper 
limit for the planet mass proposed by IQuillenl l|200qV In 
contrast, Epsilon Eri dani's disk does show clumpy structure 
|Poulton et al.l l2006). If particles are produced outside res- 
onance (as explored here) then we would infer that a mod- 
erate mass planet and low dispersion disk might account for 
the morphology. 



5 SUMMARY 

We have simulated the effect of collisions on a particle disk 
by introducing velocity perturbations into numerical inte- 
grations of particles under the influence of gravity from a 
star and a single planet. We find that the planet's chaotic 
zone boundary is an approximate location for the disk edge 
for disks that have optical depths in the range of observed 
debris disks. We find that density profiles do depend on the 
collision timescale but are more strongly dependent on the 
velocity perturbation size or velocity dispersion in the disk. 
The simulated disk morphology is axisymmetric and smooth 
unless the planet mass ratio > 10 -5 and the velocity per- 
turbation size < 0.02. For higher mass ratio planets and 
lower velocity perturbation simulations structure is seen in 
the vicinity of strong mean motion resonances at locations 
where particles are more quickly removed from the simula- 
tion due to interactions with the planet. 

From the simulations we have measured a density decre- 
ment describing the log of the density ratio of that close to 
the planet and that past the disk edge. The density decre- 
ment depends on powers of the planet mass ratio, collision 
timescale and velocity perturbation size. With a diffusion 
approximation this decrement represents a scale length that 
might in future be predicted with better understanding of 
particle dynamics near the planet. 

We have used our numerically measured density decre- 
ment to predict the mass of a planet required to truncate a 
particle disk (equation I21|l . The minimum planet mass re- 
quired depends on the square of the disk velocity dispersion, 
similar to gap opening in an accretion disk. However we find 
that a lower mass planet can truncate a low optical depth 
disk than an accretion disk with an equivalent viscosity. The 
minimum planet mass also depends on the collision rate. As 
the collision rate is related to the disk optical depth and the 
velocity dispersion related to the disk thickness, equation l21l 
can be used to estimate the mass of an unseen planet from 
observations of debris disks. 

Support for this work was in part provided by Na- 
tional Science Foundation grant AST-0406823, and the Na- 
tional Aeronautics and Space Administration under Grant 
No. NNG04GM12G issued through the Origins of Solar Sys- 
tems Program, and HST-AR-10972 to the Space Telescope 
Science Institute. 
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Figure 5a. The log of the density decrement, A, (related to the 
ratio of the density past the chaotic zone to that near the planet, 
sec equation 1 1 6 6 versus the log of the velocity perturbation size. 
Only the simulations with N c = 10 -2 collisions per orbit are 
plotted. Squares, circles and points represent simulations with 
planet mass ratio fi = 10~ 4 ,10~ 5 and 10~ 6 respectively. Also 
plotted is a linear function described by equation 1171 
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Figure 6. Disks can be truncated by planets with a mass ratio 
that lies above the lines drawn. The rr-axis is the disk velocity 
dispersion. The solid and dashed line refer to particle disks with 
a limit set by equation 1171 for N c = 10 -2 and 10 — 3 collisions 
per particle per orbit, respectively. The thin dotted line is the 
gap opening criteria expected for an a accretion disk with a = 
0.001. The dot dashed line is a limit set by requiring that the 
disk thickness be lower than chaotic zone width. A circle has been 
placed at the estimated planet mass lower limit for a planet inside 
Fomalhaut's disk edge using N c = 3 X 10 ~ 2 expected from the 
normal optical depth, t„ = 1.6 X 10~ 3 ijMarsh et alj|2005h and 
velocit y dispersion estimated from the disk edge slope u/vx ~ 
0.013 iQuillenll200e3) . 
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Table 1. List of Simulations 



Name 








CI 


itr 4 


10- 2 


0.02 


C2 


10- 5 


10~ 2 


0.02 


C3 


10~ 6 


10~ 2 


0.02 


Dl 


10- 4 


10~ 3 


0.02 


D2 


10~ 5 


10~ 3 


0.02 


D3 


10~ 6 


10~ 3 


0.02 


El 


10- 4 


10- 3 


0.04 


E2 


10- 5 


10- 3 


0.04 


E3 


10~ 6 


10~ 3 


0.04 


Fl 


10- 4 


10- 3 


0.01 


F2 


10- 5 


10- 3 


0.01 


F3 


10~ 6 


10- 3 


0.01 


Ml 


itr 4 


10~ 2 


0.04 


M2 


10- 5 


10~ 2 


0.04 


M3 


10~ 6 


10~ 2 


0.04 


Nl 


10- 4 


10~ 2 


0.01 


N2 


10- 5 


10~ 2 


0.01 


N3 


10~ 6 


10- 2 


0.01 


LI 


10- 4 


10- 2 


0.08 



The parameters describing the discussed simulations are as follows: The planet mass ratio is 
fi, whereas N c is the average number of collisions suffered by each particle per planet orbit. 
The dispersion of the azimuthal velocity kick given each collision, oa,vl v K = Zinit I Vs in units 
of the velocity of a particle in a circular orbit. The planet was in a circular orbit. During each 
collision event the radial velocity component is set to zero (e = 0). Particles were removed 
from the simulation if semi-major axis a > 1.5 or a < 0.7, eccentricity e > 0.5 or the particle 
passed within 0.1 Hill radii of the planet. Particles were generated and regenerated when 
removed with a semi-major between 1.4 and 1.5 and a mean eccentricity of ei„it- Simulations 
were run 10 s planetary orbits excepting the F series that were run twice as long and the 
LI simulation that was run half as long. These times are long because we allowed particles 
to diffuse toward the planet from a moderately distant location and we ran each simulation 
sufficiently long to achieve a steady state. 













N c = 10" 2 










□ 

o 

A 


|a=10 6 








- 





























"■■"--a... 















-2.2 



-1.8 -1.6 -1.4 

lo 9lo( e init) 



-1.2 



-1 





□ 








N c = 10" 3 






□ 




□ 

o 

A 


|a=10 6 
































A 







-2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 

lo 9l0( e init) 



